function [ lhs rhs lhsData lhsSmooth ] = HornSchunkEquations(I1, I2, x, params)
    [ Ix Iy It ] = HSGradients(I1, I2, params);
    
    Ix = Ix(:);
    Iy = Iy(:);
    It = It(:);
    
    lhsData = DataTermMatrix(Ix, Iy);
    lhsSmooth = SmoothnessTermMatrix(size(I2));
    
    lhs = lhsData - params.alpha * lhsSmooth;
    rhs = [ -Ix.*It; -Iy.*It ] + params.alpha * lhsSmooth * x;
end

function dataTerm = DataTermMatrix(Ix, Iy)    
    n = numel(Ix);
    sIx2   = spdiags(Ix.^2, 0, n, n);
    sIy2   = spdiags(Iy.^2, 0, n, n);
    sIxy   = spdiags(Ix.*Iy, 0, n, n);
    dataTerm = [ sIx2, sIxy; sIxy, sIy2 ];   
end

function sMat = SmoothnessTermMatrix(imSize)
    m = imSize(1);
    n = imSize(2);   
        
    mask = [ 0 1 0 ; 1 0 1 ; 0 1 0 ];
    [ maskRow maskCol ] = find(mask);
    
    numEntries = (4*(n-2)*(m-2)+3*2*(n-2+m-2)+4*2+n*m);
    rowInd = zeros(1, numEntries);
    colInd = zeros(1, numEntries);
    data = zeros(1, numEntries);
    
    dataCounter = 1;
    for col = 1 : n
        for row = 1 : m
            curMaskRow = maskRow + row - 2;
            curMaskCol = maskCol + col - 2;
            goodLocs = find((curMaskRow > 0) .* (curMaskCol > 0) .* (curMaskRow <=m) .* (curMaskCol <= n));
            neigInd = sub2ind(imSize, curMaskRow(goodLocs), curMaskCol(goodLocs));
            curInd = sub2ind(imSize , row, col);            
            curRowInd = repmat(curInd, 1, 1 + numel(goodLocs));            
            
            curAssignInd = dataCounter : dataCounter + numel(curRowInd) - 1;
            
            rowInd(curAssignInd) = curRowInd;
            colInd(curAssignInd) = [ curInd, neigInd' ];
            data(curAssignInd) = [-numel(goodLocs), ones(1, numel(goodLocs)) ];
            dataCounter = max(curAssignInd) + 1;
        end
    end
    
    smoothingMat = sparse(rowInd, colInd, data, m*n, m*n);
    sMat = [ smoothingMat, sparse(m*n, m*n); sparse(m*n, m*n), smoothingMat ];    
end